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Abstract 

A  finite  volume  model  of  a  solid  oxide  fuel  cell  has  been  developed.  The  model  applies  a  detailed  electrochemical  and  thermal  analysis 
to  a  tubular  SOFC  of  given  geometry,  material  properties  and  assigned  input  flows.  Electrochemical  modeling  includes  an  evaluation  of 
ohmic,  activation  and  diffusion  losses  as  well  as  a  kinetic  model  of  hydrocarbon  reactions,  based  on  most  recent  literature  experiences. 
Internal  heat  exchange  coefficients  have  been  calculated  with  a  specific  fluid-dynamic  finite  volume  analysis.  The  model  is  calibrated  on 
the  available  experimental  data  for  atmospheric  and  pressurized  tubular  SOFCs,  showing  the  capacity  of  predicting  accurately  the  SOFC 
operating  conditions.  The  model  generates  total  cell  balances  and  internal  cell  profiles  for  any  relevant  thermodynamic  or  electrochemical 
variable,  giving  the  possibility  of  discussing  the  effects  of  different  operating  conditions  on  the  internal  FC  behavior.  A  sensitivity  analysis 
is  carried  out  to  investigate  the  effects  of  different  assumptions  on  a  selection  of  key  model  parameters  involved  in  the  calculation  of  cell 
losses,  internal  heat  exchange  process  and  reforming  reactions.  Among  other  results,  it  is  shown  that  the  importance  of  the  adoption  of 
appropriate  parameters  for  the  evaluation  of  activation  polarization,  as  well  as  the  relevance  of  a  kinetic  model  for  reforming  reactions. 

©  2004  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

The  analysis  of  high  temperature  fuel  cell  power  plant 
performances  is  often  based  on  a  simplified  calculation  of 
the  cell  efficiency.  Black  box  modeling  of  the  cell  reactor, 
for  instance  with  calculation  of  cell  global  energy  balances 
and  input/output  flow  properties,  is  generally  adequate  for 
most  thermodynamic  analysis  of  complex  cycles,  e.g.  with 
integration  of  the  SOFC  with  gas  turbines  and/or  bottoming 
cycles  [1-4]. 

As  a  counterpart,  a  more  detailed  investigation  of  SOFC 
plant  performances  may  be  required  when  dealing  with 
non-standard  FC  operating  conditions  (e.g.  with  variable 
inlet  flow  conditions,  transient  or  partial  load  simulations) 
[5].  Different  approaches  are  possible  for  detailed  fuel  cells 
analysis:  finite  difference,  finite  volume  and  finite  element 
are  the  most  common  modeling  methods  adopted  in  lit¬ 
erature  [16,17].  Commercial  fluid  dynamic  computational 
codes  (CFD)  can  be  used  as  well:  in  such  cases  an  user 
defined  routine  has  normally  to  be  implemented  in  order  to 
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deal  with  electrochemical  phenomena  that  are  not  gener¬ 
ally  taken  into  account  in  commercial  CFD  codes.  The  use 
of  the  finite  volume  method  generally  yields  advantages 
in  terms  of  flexibility,  reasonable  accuracy  and  (mostly 
with  respect  to  CFD)  shorter  computational  time.  In  this 
work  we  have  defined  a  finite  volume  model  of  an  SOFC 
which  is  based  on  a  detailed  electrochemical  analysis  and 
a  fluid-dynamic  calculation  of  internal  heat  transfer  condi¬ 
tions.  Electrochemical  modeling  includes  an  evaluation  of 
ohmic,  activation  and  diffusion  losses  as  well  as  a  kinetic 
model  of  hydrocarbon  reactions,  based  on  most  recent  liter¬ 
ature  experiences.  The  model  has  been  applied  to  a  tubular 
geometry,  and  is  calibrated  on  the  available  experimental 
data  for  atmospheric  and  pressurized  tubular  SOFCs. 

The  model  calculates  cell  internal  temperature  and  flow 
composition  profiles,  fuel  and  oxidant  utilization,  cell  power 
output  and  cell  voltage  or  current  output  (depending  on  the 
calculation  option).  The  model  has  been  finally  tested  with  a 
sensitivity  analysis,  aiming  to  investigate  the  effects  of  dif¬ 
ferent  assumptions  on  a  selection  of  key  model  parameters. 

In  future  works,  the  model  will  be  applied  to  simulate 
SOFC  performances  when  integrated  in  more  complex  plant 
configurations,  enabling  to  discuss  the  effects  of  different 
operating  conditions  on  the  internal  FC  behavior. 
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Nomenclature 

A  area  (m2) 

D  mass  diffusion  coefficient  (m2/s) 

Di  j  mutual  diffusion  coefficient  of 

specie  i  in  specie  j  (m2/s) 

Dt  m  mass  diffusion  coefficient  of  specie  i  in  the 
mixture  (m2/s) 

Di  p  Knudsen  diffusion  coefficient  of  specie 
i  (m2/s) 

Lact  activation  energy  (J/mol) 

F  Faraday’s  constant:  96 487  C/mol 

h  specific  enthalpy  (J/mol) 

i  cell  current  (A) 

/o  exchange  current  (A) 

Ji  transport  rate  of  specie  i  (mol/m2  s) 

K  heat  exchange  coefficient  (W/m2  K) 

Ke q  equilibrium  constant  of  shift  reaction 

Kv  coefficient  in  methane  reforming  equation 

Lan  anode  thickness  (m) 

Lba  anodic  diffusion  path  from  bulk  flow  to 
electrode  site  (m) 

Lbc  cathode  diffusion  path  from  bulk  flow  to 
electrode  site  (m) 

Lcat  cathode  thickness  (m) 

Lei  electrolyte  thickness  (m) 

M  molecular  weight  (kg/kmol) 

n  number  of  electrons 

ni  molar  flow  of  specie  i  (mol/s) 

Nu  Nusselt  number 

p  partial  pressure 

P  pressure  (Pa) 

r  cell  radial  coordinate  (m) 

rcH4  rate  of  methane  reformation  (mol/m2  s) 
rpore  mean  radius  of  electrode  pore  (m) 

R  ohmic  resistance  (£2) 

Ro  universal  gas  constant:  8.3 14  J/mol  K 
T  temperature  (K) 

Uf  fuel  utilization 

Uox  oxidizer  utilization 

Veen  cell  voltage  (V) 

PNernst  Nernst  cell  potential  (V) 

Xt  molar  fraction  of  specie  i 

Greek  symbols 

/3  electron  transfer  coefficient 

y  pre-exponential  coefficient  in  activation 

polarization  equation 
8  thickness  (m) 

s  porosity 

Yjact  activation  polarization  (V) 

77dif  diffusion  polarization  (V) 

77 ohm  ohmic  polarization  (V) 

p  specific  resistivity  (Q  m) 

r  tortuosity 


Superscripts 

b  bulk  flow 

1  flow  at  the  electrode  site 

r  flow  at  the  reaction  site 

reac  reacted 

Subscripts 

a  air  preheating  tube 

amb  ambient 

an  anode 

c  cathodic  air 

cat  cathode 

el  electrolyte 

f  fuel 

m  gas  mixture 

s  solid  structure  (anode,  cathode,  electrolyte) 

T  preheating  tube 

tot  total 


2.  General  features 

The  model  has  been  developed  to  simulate  a  tubular  SOFC 
by  a  finite-volume  approach  [6,7].  The  fuel  cell  is  divided 
axially  in  a  user-defined  number  of  sections;  for  each  section 
the  electrochemical  and  thermal  equations  are  progressively 
solved  with  an  iterative  approach.  The  electrochemical 
model  is  first  solved  with  a  tentative  temperature  profile, 
yielding  values  of  chemical  species  fluxes,  cell  current  and 
electric  power  output.  Results  are  passed  to  the  thermal 
model,  repeating  the  process  until  convergence  is  reached 
according  to  a  user-defined  residual  error.  The  model  gen¬ 
erates  temperature  and  chemical  concentration  profiles  and 
yields  reactant  utilization,  cell  power  output  and  (depend¬ 
ing  on  the  calculation  option)  cell  voltage  or  current  output. 
Fuel  consumption  is  modeled  by  internal  reforming  and  CO 
shift  reactions  starting  with  a  fuel  composition  which  may 
include  any  combination  of  H2,  CH4,1  CO,  CO2,  H2O,  N2. 
Simulation  is  performed  with  the  following  assumptions: 

•  Stationary  conditions. 

•  Uniform  cell  voltage  along  the  cell  axis,  i.e.  equipotential 
electrode  surfaces  (Table  l).2 

•  Fluid-dynamic  conditions  in  the  reactant  channels  with 
Nusselt  number  profiles  are  discussed  in  Section  4. 1 . 

•  Heat  exchange  by  radiation  is  negligible. 

The  last  assumption  reflects  the  hypothesis  of  modeling 
a  tubular  cell  which  is  considered  to  be  bundled  to  similar 
cells  on  all  sides  inside  a  canister  of  several  hundreds  cells 
[1],  so  that  net  external  radiation  effects  may  be  neglected 


1  Higher  hydrocarbons  are  generally  cracked  in  external  prereforming 
sections. 

2  Electrode  materials  are  good  electric  conductors  (see  specific  resistivity 
values  in  Table  1). 
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Table  1 

Model  input  parameters 


Geometry  parameters 


Cell  length  (m) 

1.5  [28] 

Cell  outer  diameter  (m) 

2.2  x  10“2  [28] 

Anode  thickness  (m) 

100  x  10“6  [28] 

Cathode  thickness  (m) 

2.2  x  10“3  [28] 

Electrolyte  thickness  (m) 

40  x  10“6  [28] 

Interconnection  thickness  (m) 

85  x  10“6  [28] 

Material  properties 

Specific  resistivity  of  anode  (£2m) 

2.98  x  10“5 
exp(— 1392/TS)  [17] 

Specific  resistivity  of  cathode  (£2m) 

8.11  x  10“5  exp(600/rs) 
[17] 

Specific  resistivity  of  electrolyte 

2.94  x  10“5 

(Q  m) 

exp(10350/rs)  [17] 

Specific  resistivity  of  interconnect 

0.025  extrapolated  from 

(Q  m) 

[15,17] 

Conductivity  of  anode  (W/mK) 

2.0  [16] 

Conductivity  of  cathode  (W/mK) 

2.0  [16] 

Conductivity  of  electrolyte  (W/mK) 

2.0  [16] 

Conductivity  of  air  injection  tube 

-0.0096  x  T  +  17.892 

(W/mK) 

[22] 

Methane  reforming  (Eq.  (9)) 

Coefficient,  Kx 

8542  [9] 

Coefficient,  a 

0.85  [9] 

Coefficient,  f 

-0.35  [9] 

Ea  (kJ/mol) 

95.0  [9] 

Activation  polarization 

Activation  energy  of  anode  (kJ/mol) 

110  adapted  from  [16,19] 

Activation  energy  of  cathode  (kJ/mol) 

120  adapted  from  [16,19] 

Pre-exponential  coefficient  for  anode 

7  x  109  adapted  from 

(A/m2) 

[16,19] 

Pre-exponential  coefficient  for 

7  x  109  adapted  from 

cathode  (A/m2) 

[16,19] 

Coefficient,  m  (Eq.  (15)) 

1.0  [19] 

Diffusion  polarization 

Pore  diameter  of  anode  (m) 

^0 

H 

l 

0 

H 

X 

H 

Pore  diameter  of  cathode  (m) 

1  x  10“6  [16] 

Porosity  of  anode  (%) 

50  [16] 

Porosity  of  cathode  (%) 

50  [16] 

Tortuosity  of  anode 

3.0  [16] 

Tortuosity  of  cathode 

3.0  [16] 

thanks  to  the  geometrical  symmetry  of  the  configuration. 
At  the  cell  internal  side,  radiation  between  the  solid  cell 
structure  and  the  gas  flow  is  instead  generally  neglected  due 
to  the  dominant  effect  of  convective  heat  transfer. 

Input  data  for  the  simulation3  are: 

•  Number  of  cell  axial  sections. 

•  Cell  geometry  (length,  diameter,  electrodes  and  electrolyte 
thickness  and  electrical  properties,  injector  tube  diameter 
and  thickness). 

•  Inlet  fuel  and  air  flow  thermodynamic  properties  ( T ,  P) 
and  chemical  compositions. 

Two  calculation  options  are  then  allowed,  the  first  with 
given  cell  voltage  and  with  the  calculation  performed 


3  A  comprehensive  list  of  model  input  parameters  is  given  in  Table  1. 


straightforward,  discussed  in  the  following,  the  second  with 
given  cell  current  (or  current  density  given  the  cell  active 
area)  and  an  iterative  procedure  based  on  a  first  trial  voltage 
value. 


3.  Electrochemical  model 


The  electrochemical  model  calculates  for  each  cell  sec¬ 
tion  (Fig.  1)  the  current  power  output  and  the  molar  compo¬ 
sitions  of  cathode  and  anode  flows.  Thanks  to  the  cylindrical 
symmetry  of  the  tubular  cell  configurations,  each  section  is 
made  of  three  portions  or  finite  volumes: 

•  Fuel  (Fig.  1  shows  only  H2  for  clarity),  on  the  anode  side. 

•  Solid,  made  by  the  anode-electrolyte-cathode  structure. 

•  Oxidizer  (air),  modeled  as  an  O2/N2  mixture. 

In  each  section,  calculation  is  based  on  inlet  flow  compo¬ 
sitions  known  by  the  previous  section,  with  the  exception  of 
the  first  one  that  uses  the  assigned  input  conditions.  In  the 
anode  channel,  reforming  and  shift  reactions  are  calculated 
as  discussed  below  (Section  3.1).  Once  the  reactant  compo¬ 
sitions  are  known,  cell  current  is  calculated  by 

Tcell  =  fNernst  Oohm  —  ^?act  Odif  —  f(0  (1) 


where  Vceii  is  the  cell  voltage  and  the  Nernst  potential  VNernst 
is  calculated  by 

W,  =  £»  +  ^.n(^)+0.5l„(^-)  ,2) 

with  E°  =  1.2723  —  2.7645  x  10-4TS  the  ideal  voltage  for 
hydrogen  oxidization  at  ambient  pressure,  as  a  function  of 
temperature  at  cell  reaction  sites  [15]. 

The  ohmic,  activation  and  concentration  polarization  are 


7?  ohm  —  ^ohnfi 

(3) 

=  C°de(0  +  >?a"(0 

(4) 

»?dif  =  ^Ode(0  +  ^“(0 

(5) 

calculated  as  discussed  below  (Sections  3. 2-3. 4). 

The  non-linear  system  of  Eqs.  (l)-(5)  can  be  solved  by 
the  Muller  method  [8],  yielding  the  current  output  of  the 
considered  cell  section. 

The  number  of  moles  of  hydrogen  and  oxygen  consumed, 
as  well  as  the  generated  water  are  then 


reac  _  „reac  _ 

h2  2F  ’  °2  4  f’ 


prod 

”h2o  - 


i 

IF 


It  is  finally  possible  to  find  the  composition  of  the  flows 
given  to  the  next  cell  section  as 


nu2  ~ 
nitl  = 


h2o 


1  reac 

"h2  -  "H2  > 

J  4-  nvmd 
nU20  +  nH20 


J+1  —  J  _  wreac 
no2  ~  no2  no2  ■ 
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Fig.  1.  Electrochemical  model  principles. 


3.1.  Reforming  and  CO-shift  reactions 

Reforming  reactions  generate  Fl2  and  CO  from  the 
methane  and  steam  contained  in  the  fuel  flow: 

CH4  +  H20  =  3H2  +  CO  (8) 

Although  several  models  assume  a  simplified  approach  of 
thermodynamic  equilibrium,  a  more  correct  model  of  the 
reforming  reaction  has  to  be  based  on  a  kinetic  approach. 
Specific  literature  on  steam  reforming  in  presence  of  typi¬ 
cal  Ni/YSZ  SOFC  materials  [9-13]  suggests  the  use  of  the 
following  equation  for  the  calculation  of  the  molar  flow  rate 
of  reacted  methane  (mol/m2  s): 

rc h4  =  KrPcn4  Ph2o  exP  (  )  (9) 


where  we  have  adopted  the  values  of  Table  1  according  to 
Ref.  [9]. 

Based  on  the  section  active  surface  area  we  calculate  the 
number  of  reacted  moles  of  methane  (Veac)  and  the  fuel 
components  mole  number  after  reforming  (Anew)  as 


new  ,  old  reac 

nCK4  ~  nCK4  nCH4’ 


n 


new 

H2 


H-  3  n 


reac 

CH4 


n 


new 

CO 


X7reac 

nCO  ^  ^CH4’ 


(10) 


based  on  the  original  fuel  composition  (ftold). 

Final  fuel  composition  is  modeled  by  the  shift  reaction 


CO  +  H20  =  H2  +  C02 


(11) 


Which  is  considered  to  reach  local  equilibrium  [14-16] 
with  an  equilibrium  constant  depending  on  temperature 


Cathode 


Fig.  2.  (a)  Cell  current  path;  (b)  cell  equivalent  electric  circuit. 
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^CQ2^H2 

^co^h2o 


=  exp 


3.961 


3.2.  Cell  ohmic  resistance 


(12) 


Calculation  of  the  cell  ohmic  resistance  (Eq.  (3))  is  based 
on  the  cell  electric  model  of  Fig.  2a  and  b  where  current 
flows  across  interconnections,  anode,  electrolyte  and  cathode 
under  the  cell  circumferential  potential  difference  [17].  Cell 
current  path  is  simulated  with  the  hypothesis  of  uniform 
axial  voltage  (equipotential  electrodes).4 

The  cell  equivalent  ohmic  resistance  [14,17]  depends  on 
the  anode,  cathode  and  electrolyte  resistances  RCSLt ,  Re\ 

which  are  calculated  according  to  the  second  Ohm’s  law: 


(13) 


where  Rf  is  the  ohmic  resistance  in  each  share  of  the  equiv¬ 
alent  circuit  of  Fig.  2b,  A;  the  respective  area  of  the  sec¬ 
tion  where  the  current  flows,  5;  the  corresponding  current 
flow  length  and  pi  is  the  corresponding  material  resistivity, 
calculated  with  temperature-dependent  relations  (Table  1). 
Thanks  to  the  cell  symmetry,  the  total  ohmic  resistance  is 
calculated  by  two  identical  resistance  paths  in  parallel  ar¬ 
rangement  [18]. 


3.3.  Activation  polarization 


The  development  of  electrochemical  reactions  requires 
overcoming  an  activation  energy  barrier,  yielding  a  po¬ 
tential  loss  defined  as  activation  polarization  with  the 
Butler- Volmer  relations  [16,17,19,20]: 


i  =  i  o 


exp  (  (3 


nFr ] 


act 


RgTs 


exp 


-(l-j8) 


nFr) 


act 


j 


(14) 


where  j3  is  the  electronic  transfer  coefficient  and  z'o  the  ex¬ 
change  current  density  that  can  be  calculated  as 


^0,an  —  Kan 


PU2  \  { ^H2o\ 


m 


/Aimb 


Pamb 


exp 


D  act,  an 
RgTs 


A 


(15) 


/  P o2 

^0,cat  —  Xcat  I 

\  i^amb 


act,  cat 

RgTs 


A 


(16) 


Values  for  yan,  ycat,  Fact,  an,  Fact, cat,  m  have  been  found  in 
literature  [16,19];  however,  due  to  the  particularly  wide  dis¬ 
persion  of  values  found  for  these  parameters,  it  is  important 
to  consider  the  effect  of  different  assumptions  on  their  val¬ 
ues,  as  discussed  in  the  last  section  of  this  work  (sensitivity 
analysis). 


4  It  should  be  noted  that,  despite  the  circumferential  current  path,  current 
density  for  tubular  SOFCs  is  conventionally  expressed  with  reference  to 
the  cell  external  area  (cylindrical  surface). 


If  the  SOFC  works  in  low  activation  polarization  condi¬ 
tions,  Eq.  (14)  can  be  approximated  in  a  linear  form: 


7?  act  — 


Fg  Tsi 
nFio 


(17) 


If  SOFC  operation  features  high  polarization,  it  is  possible 
to  neglect  the  second  term  of  Eq.  (14)  and  write  the  relation 
known  as  Tafel’s  law: 


7?  act  — 


^Mn- 

nF/3  i0 


(18) 


where  the  value  of  j3  is  chosen  to  keep  continuity  of  Eqs.  (17) 
and  (18). 

According  to  the  approach  used  in  [16],  low  versus  high 
polarization  conditions  are  determined  by 


F?7act  1 

-  <  1 

nRgTs 


(19) 


3.4.  Diffusion  polarization 


Calculation  of  cell  reversible  potential  in  Eq.  (2)  is  based 
on  the  average  “bulk  flow”  reactant  chemical  composition. 
Due  to  reactant  consumption  at  the  cell  surface  a  diffusion 
mass  transfer  process  occurs,  where  the  real  concentration 
at  cell  reaction  sites  and  the  corresponding  cell  potential  is 
lower.  The  difference  between  “bulk”  potential  and  reaction 
site  potential  is  called  diffusion  polarization  and  is  given  by 
two  terms,  related  to  the  anode  and  cathode  sides: 


mu  = 


¥1 

IF 


In 


*H2*H2o\ 

*h2o*hJ 


RpTs 

+ 

4  F 


Vo2 


(20) 


where  bulk  and  reaction  site  concentrations  are  indicated 
with  apex  “b”  and  “r”,  respectively. 

Calculation  of  the  reactant  molar  fraction  at  cell  reaction 
sites  is  carried  out  as  discussed  in  Appendix  A. 


4.  Heat  exchange  model 

The  electrochemical  analysis  of  the  fuel  cell  carried  out 
in  Sections  3. 1-3.4  is  strongly  influenced  by  temperature. 
The  detailed  fuel  cell  simulation  then  requires  an  analysis  of 
the  internal  temperature  profile.  The  heat  exchange  model 
is  based  on  the  same  axial  sections  defined  for  the  electro¬ 
chemical  model;  each  section  is  then  divided  into  four  finite 
volumes,  as  shown  in  Fig.  3: 

•  air  preheating  inside  the  injection  tube, 

•  cathode  air, 

•  solid  structure  (anode,  cathode,  electrolyte), 

•  fuel  mixture. 

On  each  of  the  four  finite  volumes  it  is  carried  out  an 
energy  balance,  yielding  4xw  equations  with  4  xn  unknown 
temperatures  for  n  cell  sections;  temperature  is  considered 


118 


S.  Campanari,  P.  lorn/  Journal  of  Power  Sources  132  (2004)  113-126 


Fig.  3.  Heat  exchange  model  configuration. 


uniform  inside  each  volume  [7].  Energy  balances  for  the 
generic  ith  section  are  discussed  below. 


(a)  Air  preheating  volume 


n 


o2  no2,Ta  +  n 

-  n02hb2, Ta 


N>N2!Ta  +  KtAT(TI 

~  «N2^N2,Ta  =  0 


(21) 


Temperature  is  influenced  by  heat  exchange  with  cath¬ 
ode  airflow  according  to  a  global  heat  exchange  coeffi¬ 
cient  Kj,  which  includes  convective  heat  transfer  inside 
and  outside  the  injection  tube  and  conductive  heat  trans¬ 
fer  across  the  tube  thickness  [21].  Thermal  conductivity 
of  the  alumina  tube  is  expressed  as  a  function  of  tem¬ 
perature  according  to  [22]. 

(b)  Cathode  air  volume 


Heat  is  exchanged  with  cathode  air,  with  the  fuel  mix¬ 
ture  and  with  the  two  adjoining  solid  volumes,  T^sf  is  the 
convective  heat  exchange  coefficient  between  the  sold 
structure  and  the  fuel  mixture;  K$  accounts  for  conduc¬ 
tive  heat  transfer  between  two  adjoining  solid  volumes, 
calculated  based  on  anode,  cathode  and  electrolyte  av¬ 
erage  thermal  conductivity  [16]. 

(d)  Fuel  volume 


n 


i— 1 


nX,Tf 


+  n 


i 

H20,reac 


h 


H20,Tf 


x=h2,h2o,ch4,co,co2,n2 


+  KSFAsF(Ti  -  t‘)  - 

x=h2,h2o,ch4,co,co2,n2 


x  nlxhlx  Tf 


nH2,reac^H2,Tf  “  ^ 


(24) 


+  nN2l  ^n21,Tc  +  At  Ax  (7^  - 
+  ASCAsc (^  -  Tlc)  -  ^b2,reac^02,Tc 
~^N2^N2,Tc  =  ^ 


Tq) 

-  nb2*02,Tc 

(22) 


Temperature  is  influenced  by  heat  exchange  with  the  air 
preheating  volume  and  with  the  solid  structure  volume. 
Ksc  is  the  convective  heat  transfer  coefficient  between 
cathode  air  and  the  solid  structure. 

(c)  Solid  structure  volume 


nb2,reac^02,Tc  +  nH2,reac^H2,Tf  +  AscAsc^  ^s) 

+  KSfASf(T)  ~  T‘)  +  KSAS(T‘+1  -  T‘) 

+  KsAsiTp1  ~  Tj)  -  n^_0  pmdh'H20  Ts  -  Wel  =  0 

(23) 


Heat  is  exchanged  with  the  solid  structure  with  the  same 
approach  used  in  the  other  volumes. 

Molar  flows  of  all  chemical  species  in  Eqs.  (21)-(24) 
are  known  by  the  electrochemical  model;  specific  en¬ 
thalpies  are  expressed  as  a  linear  function  of  tempera¬ 
ture: 


hi  =  a  +  bT  (25) 

with  coefficients  a,  b  given  by  interpolation  of  available 
enthalpy  data  [23]. 

Heat  exchange  coefficients  are  calculated  as  a  function  of 
gas  and  material  thermal  conductivities  (expressed  by  third 
grade  polynomial  functions  of  temperature  [24]),  depending 
on  the  flow  conditions. 
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Fig.  4.  Velocity  field  (m/s)  in  the  region  inside  the  cell  closed  end. 


tions  along  the  majority  of  the  fuel  cell.  They  indicate  that 
the  fluid  dynamic  field  respects  laminar  flow  conditions  at 
the  anode  side  (external  tube  wall,  Re  =  50),  while  the  air 
flow  at  cathode  side  (internal  tube  wall)  and  inside  the  in¬ 
jection  tube  is  relatively  closer  to  a  transition  to  turbulent 
flow  (Re  =  75 0-2700). 5  For  the  three  channels  the  result¬ 
ing  average  Nusselt  number  adopted  for  the  solution  of  the 
heat  exchange  model  are  then  equal  to  4.2  (anode  side),  5.5 
(cathode  side)  and  1 1  (injection  tube).  Nusselt  numbers  have 
been  calculated  from  total  heat  exchange  coefficients,  with 
reference  to  wall  temperatures  and  channels  adiabatic  mix¬ 
ing  temperature. 


4.1.  Gas  flow  fluid- dynamic  conditions 

High  temperature  fuel  cell  models,  which  can  be  found  in 
the  available  literature,  are  generally  based  on  the  assump¬ 
tion  of  a  constant  Nusselt  number  for  the  calculation  of  cell 
internal  heat  exchange  process.  A  Nusselt  number  equal  to  4, 
typical  of  low  speed  laminar  flow  conditions  is  for  instance 
assumed  both  for  fuel  and  air  side  heat  exchange  [16,19]. 

Aiming  to  better  focus  this  subject,  we  have  carried  out  a 
detailed  investigation  of  the  gas  flow  conditions  inside  and 
outside  the  fuel  cell  (air  and  fuel  flows)  with  a  finite  volume 
analysis  based  on  a  computational  fluid  dynamics  (CFD) 
software  (Fluent®). 

The  CFD  simulation  has  been  performed  with  the  same 
air  and  fuel  flow  inlet  conditions  used  for  the  model  calibra¬ 
tion  of  Section  5.  A  uniform  heat  generation  is  considered 
into  the  solid  structure  in  order  to  simulate  the  exothermic 
electrochemical  reactions;  a  symmetry  condition  is  applied 
to  the  external  boundary,  assumed  at  a  radial  distance  of 
16  mm  from  the  cell  axis  [14].  The  simulation  code  solves 
fluid  velocity  and  temperature  fields  with  a  segregated  up¬ 
wind  second-order  method;  following  several  simulations, 
the  k-co  viscosity  model  has  been  used. 

The  considered  cell  geometry  is  indicated  in  Table  1, 
while  Fig.  4  shows  a  detailed  example  of  the  velocity  field 
in  the  region  of  the  cell  closed  end,  where  the  model  grid 
needs  particular  refinement. 

Results  for  the  Reynolds  number  are  shown  in  Fig.  5 
for  a  cell  radial  section  representative  of  the  flow  condi- 


4.2.  Boundary  conditions 

As  the  cell  is  divided  into  n  sections,  boundary  conditions 
are  applied  to  the  external  limits  of  cells  i  =  0  and  i  =  n  + 1 . 

The  flow  transport  process  of  the  air  and  fuel  flow  consti¬ 
tutes,  mathematically  speaking,  a  parabolic  problem,  where 
the  solution  for  each  section  is  not  influenced  by  downward 
sections;  its  solution  then  requires  an  initial  condition  for 
the  inlet  flows  (Ta?in,  7f?in) 

Tfln  +  1)  =  Ta,in  (26) 

7f  (0)  =  7f,in  (27) 

Heat  exchange  process  inside  the  solid  structure  has  instead 
a  mathematical  elliptic  nature,  and  the  solution  for  each 
section  depends  on  both  downward  and  upward  solutions. 
Two  boundary  conditions  for  both  cell  extremity  are  then 
required:  it  has  been  considered  here  that:  (i)  first  solid  cell 
temperature  can  be  expressed  as  an  average  between  fuel  and 
air  temperature,  weighted  according  to  the  results  of  the  heat 
transfer  analysis  discussed  in  Section  4.1  and  keeping  into 
account  the  endothermic  effects  of  the  reforming  reactions, 
and  (ii)  heat  exchange  is  zero  at  the  outlet  extremity.  The 
conditions  are  then  expressed  as 

rs(0)  =  7f(0)  X  0.4  +  rc(0)  X  0.6  (28) 

Ts(n  +  l)  =  Ts(n)  (29) 

A  further  condition  is  applied  to  express  continuity  between 
the  airflow  exiting  the  injection  tube  and  entering  the  cath¬ 
ode: 
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Fig.  5.  Reynolds  number  in  a  central  section  of  the  fuel  cell. 


Tc(  0)  =  ra(0)  (30) 

4.3.  Solution  algorithm 

The  electrochemical  model  allows  being  calculated  one 
section  at  a  time,  as  the  solution  for  each  section  gives  in¬ 
put  conditions  for  the  following.  As  already  mentioned,  the 
thermal  model  instead  generates  4  xn  equations  with  4  xn 
unknown  temperatures.  With  an  upwind  solution  approach 


5  Calculation  of  Reynolds  number  is  based  on  hydraulic  diameter  for 
the  anode  and  cathode  channels. 
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Table  2 

Data  for  model  calibration 


Plant  A 

Plant  B 

SOFC  module  power  output  (kWdc) 

120.7 

267.5 

Cell  number 

1152 

1704 

Operating  pressure  (bar) 

1.05 

3.5 

Fuel  mass  flow  (kg/s) 

0.03578 

0.08232 

Air  mass  flow  (kg/s) 

0.35 

0.634 

Cell  voltage  (V) 

0.69 

0.639 

Cell  inlet  fuel  temperature  (°C) 

550 

587 

Cell  inlet  air  temperature  (°C) 

831 

775 

Fuel  composition  at  cell  inlet 

(molar  fraction) 

h2 

0.258 

0.226 

h2o 

0.284 

0.334 

ch4 

0.11 

0.131 

CO 

0.057 

0.057 

co2 

0.228 

0.241 

n2 

0.063 

0.011 

[6]  the  model  generates  a  linear  three-diagonal  system  of 
equations,  which  can  be  solved  by  the  Gauss-Seidel  method 
without  heavy  computational  requirements. 

For  a  number  of  cell  sections  equal  to  100  (a  number 
far  sufficient  to  correctly  model  the  cell  behavior  as  shown 
in  Table  9)  a  total  calculation  time  in  the  range  of  few 
10  s  is  required  to  reach  convergence  on  an  average  power 
PC. 


5.  Model  calibration 

A  model  calibration  has  been  performed  utilizing  avail¬ 
able  data  for  two  tubular  SOFC  modules:  (i)  an  SOFC 
module  operated  in  a  lOOkW  atmospheric  prototype  plant 
[25-30]  and  (ii)  an  SOFC  module  manufactured  for  the  in¬ 
tegration  with  a  recuperated  gas  turbine  cycle  in  a  300  kW 
“hybrid  plant”,  where  the  SOFC  is  operated  under  pressur¬ 
ized  conditions  [31-33]. 

The  input  calculation  parameters  are  shown  in  Table  1, 
with  reference  to  the  model  discussion  given  above  and  in 
Appendix  A.  Both  modules  are  based  on  the  same  tubular 
SOFC  units,  so  that  the  model  input  parameters  are  cali¬ 
brated  with  the  same  values  for  both  cases. 


Fig.  6.  Temperature  profile.  In  all  cases  on  the  horizontal  axis  it  is  shown 
the  cell  axial  coordinate  starting  from  the  reactant  inlet  side  (tube  closed 
end,  see  Fig.  1). 

Calibration  data  for  both  power  plants  are  shown  in 
Table  2;  first  data  are  based  on  a  literature  research  while 
flow  temperatures  and  composition  are  based  on  previous 
modeling  works  [4,25]. 

Results  of  calibration  are  shown  in  Table  3.  They  show 
that  the  model  input  parameter  set  in  Table  1  allows  to  cor¬ 
rectly  match  the  two  SOFC  performances;  the  agreement  be¬ 
tween  calculated  and  expected  performances  is  particularly 
good  for  the  first  SOFC,  while  keeping  reasonable  accuracy 
also  for  the  second  case  (maximum  error  close  to  3%). 

Better  accuracy  could  be  reached  by  further  custom¬ 
tailoring  some  of  the  most  effective  and  uncertain  model 
input  parameters  (see  also  the  sensitivity  analysis  in  last 
section  of  this  work),  for  instance  those  involved  by  ac¬ 
tivation  losses  (Eqs.  (15)  and  (16)),  as  well  as  by  having 
more  accurate  data  for  the  inlet  flow  temperatures  and 
compositions,  which  are  generally  kept  proprietary. 

6.  Model  results 

The  model  results  are  primarily  constituted  by  total  cell 
balances  and  by  cell  internal  profiles  for  any  relevant  ther¬ 
modynamic  or  electrochemical  variable.  The  first  results  are 
anticipated  by  Table  3,  and  they  include  cell  power  output, 
reactant  utilization  and  cell  voltage  or  total  current  (depend¬ 
ing  on  the  calculation  option).  Cell  internal  profiles  include 
the  examples  of  Figs.  6-10.  In  all  cases  on  the  horizontal 


Table  3 

Calibration  results 


Plant  A 

Plant  B 

Calculated 

Expected 

Error  (%) 

Calculated 

Expected 

Error  (%) 

Single  cell  power  (W) 

103.1 

104.8 

1.6 

161.7 

157.0 

3.0 

Current  density  (A/m2) 

1792 

1800 

0.4 

3034 

3000 

1.1 

Fuel  utilization3  (%) 

68.6 

69.0 

0.6 

70.9 

69.0 

2.9 

Air  utilization  (%) 

17.5 

17.8 

1.7 

24.2 

23.8 

1.7 

Current  density  is  calculated  with  a  cell  active  area  equal  to  834  cm2  [29]. 
a  Single  passage,  without  fuel  recycling. 
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Fig.  7.  Flow  chemical  composition  profiles. 

axis  it  is  shown  the  cell  axial  coordinate  starting  from  the 
reactant  inlet  side  (tube  closed  end,  see  Fig.  1). 

Fig.  6  shows  the  temperature  profile  for  the  fuel  flow  (7f ), 
the  solid  wall  (Ts),  the  cathode  air  flow  (rc)  and  the  injec¬ 
tion  tube  internal  air  flow  (Ta,  dashed  line).  The  fuel  flow, 
as  well  as  the  solid  wall,  is  progressively  heated  up  by  the 
cell  heat  generation,  but  the  temperature  rise  is  slower  in  the 
first  portion  of  the  fuel  flow  path  due  to  the  hydrocarbon  re¬ 
forming  reactions.  In  this  area,  the  injected  air  and  the  cath¬ 
ode  airflow  both  act  as  thermal  energy  storage,  preventing 
the  formation  of  “cold  spots”  or  any  excessive  decrease  of 
temperature. 

In  the  region  close  to  the  cell  inlet,  the  solid  structure  is 
heated  by  the  hot  cathode  air  flow  and  cooled  by  the  fuel 
flow.  Fuel  flow  temperatures  do  not  allow  in  principle  a 
fast  reforming  reaction,  so  that  hydrogen  consumption  by 
electrochemical  reactions  causes  a  local  minimum  in  H2 
concentration  (Fig.  7)  in  the  area  at  about  5  cm  from  the  cell 
inlet.  As  soon  as  the  fuel  flow  is  heated  up,  the  prevailing 
phenomenon  becomes  the  reforming  reaction  leading  to  an 
increase  in  H2  concentration  up  to  a  maximum  at  about 
40  cm  from  the  cell  inlet. 

The  progression  of  reforming  reactions  is  shown  also 
by  the  methane  concentration,  with  CH4  disappearing  after 
about  one  third  of  the  cell  length.  Hydrogen  and  CO  con¬ 
sumption  leaves  the  way  to  a  rapid  formation  of  oxidization 
products  (CO2,  H2O)  when  the  fuel  cell  oxidization  reac¬ 
tions  become  dominant. 


Fig.  9.  Ohmic  and  activation  polarization  profiles. 


At  the  cell  outlet  (right  side  in  Fig.  6),  the  fuel  and  cath¬ 
ode  flow,  as  well  as  the  solid  structure,  show  a  temperature 
decrease  due  to  the  heat  transfer  to  the  incoming  airflow 
(dashed  line)  which  is  progressively  heated  up  inside  the  in¬ 
jection  tube. 

Fig.  8  shows  a  progressive  decrease  of  current  output  and 
of  cell  reversible  potential  (Nernst  potential)  after  the  first 
half  of  fuel  cell:  following  the  end  of  the  reforming  reac¬ 
tions,  the  hydrogen  moles  consumed  by  the  fuel  cell  reac¬ 
tions  are  no  longer  balanced  by  the  hydrogen  moles  gen¬ 
erated  by  reforming,  while  the  steam  fraction  continuously 
increases.  The  result  is  a  progressive  decrease  of  hydrogen 
concentration,  yielding  a  Nernst  potential  and  current  output 
decrease.  Moreover,  low  value  of  current  density  in  the  first 
half  are  caused  by  solid  structure  low  temperature  (Fig.  6) 
yielding  high  activation  and  ohmic  polarization  losses. 

The  three  kinds  of  polarization  losses,  which  contribute 
to  the  decrease  of  the  cell  potential  from  the  ideal  Nernst 
value  to  the  real  value,  are  shown  in  Figs.  9  and  10.  Ac¬ 
tivation  losses,  both  at  the  anode  and  cathode  side,  play  a 
significant  role  in  the  first  portion  of  the  fuel  cell,  where  the 
average  temperatures  are  relatively  low.  Ohmic  losses  are 
proportional  to  the  cell  current  density  shown  in  Fig.  8  and 
are  as  well  influenced  by  lower  temperature  in  the  first  por¬ 
tion  as  predicted  by  the  assumption  made  on  cell  resistivity 
(Table  1).  Diffusion  losses  (Fig.  10)  are  negligible  with  re¬ 
spect  to  the  other  losses. 


E 

5 


0.08 


Fig.  8.  Current  density  and  reversible  cell  voltage  profiles. 


Fig.  10.  Diffusion  polarization  profiles. 
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Table  4 


Activation  polarization  sensitivity  analysis  (reference  case  as  in  Table  1) 


Em  (kJ/kg) 

90 

100 

105 

110 

112 

115 

118 

Cell  power  (W) 

107.7 

106.8 

105.7 

103.1 

100.8 

75.8 

41.5 

i  cell  (A/m2) 

1871 

1856 

1837 

1792 

1750 

1316 

721 

Uf  (%) 

71.7 

71.1 

70.3 

68.6 

67.0 

50.4 

27.6 

uox  (%) 

18.1 

18.1 

18.0 

17.5 

17.1 

12.9 

7.1 

Refend  (mm) 

420 

435 

465 

525 

585 

1065 

1455 

V  an 

3  x  109 

4  x  109 

5  x  109 

7  x  109 

o 

O 

X 

H 

5  x  1010 

1  X  1011 

Cell  power  (W) 

52.5 

88.8 

99.0 

103.1 

105.1 

107.6 

107.9 

i  cell  (A/m2) 

906 

1542.5 

1719.6 

1792 

1826 

1870 

1874.6 

Uf  (%) 

34.7 

59.1 

65.8 

68.6 

69.9 

71.6 

71.7 

Uox  (%) 

8.9 

15.1 

16.8 

17.5 

17.8 

18.3 

18.3 

Refend  (mm) 

1350 

840 

630 

525 

480 

420 

405 

£cat  (kJ/kg) 

105 

110 

115 

120 

125 

128 

130 

Cell  power  (W) 

105.2 

104.9 

104.4 

103.1 

99.6 

79.0 

49.3 

i  cell  (A/m2) 

1828 

1823 

1813 

1792 

1731 

1373 

857 

Uf  (%) 

70.0 

69.8 

69.4 

68.6 

66.3 

52.6 

32.8 

Uox  (%) 

17.9 

17.8 

17.7 

17.5 

16.9 

13.4 

8.4 

Refend  (mm) 

465 

480 

495 

525 

615 

1020 

1395 

Yc  at 

2  x  109 

3  x  109 

5  x  109 

7  x  109 

o 

O 

X 

H 

5  x  1010 

1  X  1011 

Cell  power  (W) 

52.7 

94.5 

101.7 

103.1 

104.0 

105.3 

104.5 

i  cell  (A/m2) 

915 

1642 

1766 

1792 

1808 

1830 

1833 

Uf  (%) 

35.0 

62.9 

67.6 

68.6 

69.2 

70.1 

70.2 

Uox  (%) 

8.9 

16.0 

17.3 

17.5 

17.7 

17.9 

17.9 

Refend  (mm) 

1350 

735 

570 

525 

510 

465 

465 

Table  5 

Ohmic  polarization  sensitivity  analysis  (reference  case  as 

in  Table  1) 

Pint  (Qm) 

0.005 

0.01 

0.02 

0.025 

0.03 

0.05 

0.1 

Cell  power  (W) 

113.6 

111.3 

106.1 

103.1 

99.9 

84.8 

56.4 

i  cell  (A/m2) 

1974 

1935 

1834 

1792 

1736 

1474 

980 

Uf  (%) 

75.6 

74.1 

70.6 

68.6 

66.5 

56.4 

37.5 

Uox  (%) 

19.3 

18.9 

18.0 

17.5 

17.0 

14.4 

9.6 

Refend  (mm) 

360 

390 

480 

525 

585 

825 

1215 

7.  Sensitivity  analysis 

A  sensitivity  analysis  on  the  results  for  the  first  SOFC 
module  is  proposed,  aiming  to  highlight  the  effect  of  differ¬ 
ent  polarization  losses  on  cell  performance. 

Activation  and  ohmic  polarization  are  considered  in 
Tables  4  and  5  where  the  input  data  for  the  default  case  are 
given  in  Table  1. 

Anode  and  cathode  activation  polarization  sensitivity 
analysis  is  performed  by  varying  four  operating  parameters 
(Eqs.  (15)  and  (16)):  activation  overpotential  pre-exponential 
coefficient  and  activation  energy  of  anode  and  cathode  yan, 
Xcat,  £act,an,  ^act,cat-  These  values  are  varied  in  the  relatively 
wide  range  that  can  be  found  in  literature:  activation  energy 
ranges  between  100-140  kJ/kg  for  anode  and  117-160  kJ/kg 
for  cathode,  while  for  the  pre-exponential  values  in  the  or¬ 
der  of  magnitude  from  108  to  1010  A/m2  have  been  found 
[16,19].  All  the  other  input  data  are  assumed  as  in  Table  1. 

Increasing  anode  and  cathode  activation  energy  of  about 
5-10  kJ/kg  leads  to  a  progressively  steeper  decrease  of  cell 
power  and  current  output,  as  well  as  reactant  utilization  fac¬ 
tors.  This  is  due  to  the  exponential  nature  of  the  equations 


where  activation  energy  is  the  exponent  argument.  The  cell 
axial  position  where  the  reforming  reactions  are  completed 
(Refend)  becomes  more  distant  from  the  cell  inlet.  The  cell 
performance  is  completely  upset  for  an  activation  energy  in¬ 
crease  in  the  range  of  15-20  kJ/kg. 

A  similar  situation  arises  for  a  two  orders  of  magnitude 
decrease  of  activation  overpotential  pre-exponential  coeffi¬ 
cient,  from  1011  to  109  A/m2.  These  results  suggest  the  ex¬ 
treme  importance  of  a  correct  calibration  of  these  parame¬ 
ters  for  an  accurate  fuel  cell  modeling.  As  shown  by  Table  6, 


Table  6 

Overall  effect  of  diffusion  and  activation  losses  (reference  case  as  in 
Table  1) 


Reference 

Zero  diffusion 

losses 

Zero  activation 
polarization 

Cell  power  (W) 

103.166 

103.194 

109.5 

i  cell  (A/m2) 

1792.76 

1793.24 

1901.9 

Uf  (%) 

68.634 

68.652 

72.80 

uox  (%) 

17.518 

17.522 

18.61 

Refend  (mm) 

525 

525 

360 
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Table  7 

Overall  effect  of  Nusselt  number  (reference  case  as  in  Table  1) 


Reference 

Nu  =4 

Nu  =1 

Nu  =  10 

Nu  =  13 

Nu  15 

Cell  power  (W) 

103.1 

104.8 

103.3 

100.8 

92.0 

51.9 

i  cell  (A/m2) 

1792 

1821 

1795.9 

1751.2 

1598 

902 

Uf  (%) 

68.6 

69.7 

68.8 

67.0 

61.2 

34.5 

t/ox  (%) 

17.5 

17.8 

17.5 

17.1 

15.6 

8.8 

Refend  (mm) 

525 

510 

510 

585 

840 

1455 

Tnax  (°C),  position  (mm) 

945  at  1215 

937  at  1380 

940  at  1245 

940  at  1230 

907  at  1350 

807  at  outlet 

Tm in  (°C),  position  (mm) 

698  at  inlet 

708  at  inlet 

688  at  inlet 

661  at  inlet 

612  at  inlet 

557  at  inlet 

Table  8 

Overall  effect  of  different  methane  reforming  assumptions  (reference  case  as  in  Table  1) 

Reference 

Ref.  [16] 

Ref.  [17] 

Cell  power  (W) 

103.1 

103.0 

101.7 

i  cell  (A/m2) 

1792 

1790 

1767 

Uf  (%) 

68.6 

68.5 

67.7 

t/ox  (%) 

17.5 

17.5 

17.3 

Refend  (mm) 

525 

615 

495 

Tnax  (°C),  position  (mm) 

945  at  1215 

945  at  1215 

939  at  1260 

Tinin  (°C),  position  (mm) 

698  at  inlet 

695  at  inlet 

677  at  inlet 

Table  9 

Effect  of  different  meshing  of  the  finite  volume  model 

Number  of  cells 

10 

50 

100 

200 

500 

1000 

Cell  power  (W) 

106.5 

103.7 

103.1 

102.9 

102.8 

102.7 

i  cell  (A/m2) 

1851 

1801 

1792 

1788 

1786 

1785 

Uf  (%) 

70.9 

69.0 

68.6 

68.5 

68.4 

68.4 

t/ox  (%) 

18.1 

17.6 

17.5 

17.5 

17.5 

17.4 

Refend  (mm) 

300 

480 

525 

562 

588 

600 

Tnax  (°C),  position  (mm) 

944  at  1200 

945  at  1230 

945  at  1215 

946  at  1215 

947  at  1212 

948  at  1220 

7^  (°C),  position  (mm) 

703  at  inlet 

703  at  inlet 

698  at  inlet 

698  at  inlet 

699  at  inlet 

700  at  inlet 

rather  than  assuming  very  bad  values  for  those  parameters, 
it  could  be  better  to  simulate  the  fuel  cell  totally  neglecting 
the  activation  loss:  the  counterpart  to  this  choice  would  be 
tolerating  performance  errors  in  the  range  of  6-7%  on  power 
output  and  reactant  utilization  factors,  and  would  also  be  re¬ 
nouncing  to  correctly  model  the  internal  profile  of  reforming 
reactions.  Nevertheless,  it  should  be  noted  that  since  activa¬ 
tion  polarizations  are  strongly  temperature  dependent,  they 
are  extremely  important  in  detecting  bad  cell  operating  con¬ 
ditions  which  may  occur  at  low  cell  temperature,  and  they 
in  general  should  be  considered  in  the  evaluation  of  cell  per¬ 
formance. 

Ohmic  polarization  parametric  analysis  is  performed  by 
varying  the  interconnect  resistivity,  due  to  the  particular  wide 
range  of  values  which  are  proposed  in  literature  (for  instance 
ranging  from  0.05  Qm  [17]  to  0.01  Qm  [15]). 

The  effect  of  a  progressive  increase  of  resistivity  is  an 
almost  linear  decrease  of  power  output  and  current  density, 
as  well  as  of  fuel  and  air  utilization. 

At  the  high  operating  temperatures  typical  of  SOFCs, 
thanks  also  to  the  extreme  thinness  of  cell  active  lay¬ 
ers,  diffusion  polarization  has  a  very  small  effect  on  cell 


performance.  As  anticipated  by  the  graphical  profiles  of 
the  cell  losses,  Table  6  shows  clearly  that  calculated  val¬ 
ues  are  only  slightly  influenced  when  diffusion  losses  are 
neglected. 

Table  7  deals  with  a  sensitivity  analysis  conducted  on 
different  assumption  of  Nusselt  number.  Starting  from  the 
conditions  of  Table  1  in  the  first  column,  we  have  assumed 
a  constant  value  of  Nu  =  4-15  in  all  fuel  cell  channels 
in  the  other  columns.  It  should  be  noted  that  there  is  a 
little  difference  with  respect  to  the  reference  case  when  a 
global  laminar  flow  behavior  (Nu  =  4)  is  assumed.  Increas¬ 
ing  Nusselt  number  over  10  increases  the  heat  exchange  ef¬ 
fectiveness,  leveling  the  cell  temperature  at  lower  values, 
especially  in  the  first  portion  in  the  direction  of  the  fuel 
flow;  this  situation  significantly  slows  down  the  reform¬ 
ing  reactions  and  leads  to  a  progressive  decay  in  the  cell 
performance. 

In  Table  8  reference  case  is  compared  with  different  as¬ 
sumptions  of  the  parameters  used  for  the  simulation  of  the 
reforming  reaction  kinetics  (Eq.  (9)).  Results  in  the  third 
column  refer  to  value  of  coefficients  chosen  according  to 
[16]  where  a  first-order  dependence  on  methane  partial 
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pressure  is  assumed,  while  contribution  on  steam  partial 
pressure  is  neglected  (i.e.  value  of  a  =  1  and  /3  =  0  are 
assumed  in  Eq.  (9)).  Data  shown  in  the  last  column  refers 
to  simulation  performed  with  a  =  1  and  /3  =  —1.25  as 
proposed  in  [17]. 6  As  shown  in  Table  8,  slight  differences 
exist  with  respect  to  reference  case,  therefore  both  [16,17] 
methane  reforming  approaches  give  a  realistic  representa¬ 
tion  of  the  molar  flow  rate  of  reacted  methane.  Conversely, 
simulation  with  thermodynamic  equilibrium  approach  gives 
way  to  misleading  results,  yielding  a  too  fast  methane 
consumption  at  the  inlet  portion  of  the  cell:  the  effect 
is  an  abrupt  decrease  of  cell  temperature  which  magni¬ 
fies  the  polarization  losses  and  completely  upsets  the  cell 
simulation. 

Finally,  Table  9  shows  the  effect  of  different  meshing  on 
the  results  of  the  simulation.  It  is  evident  that  a  subdivision  of 
100-200  finite  volume  sections  leads  to  a  correct  evaluation 
of  the  cell  behavior. 


8.  Conclusions 

In  this  work  we  have  defined  a  finite  volume  model 
of  an  SOFC  which  is  based  on  a  detailed  electrochem¬ 
ical  analysis  and  a  fluid-dynamic  calculation  of  internal 
heat  transfer  conditions.  Electrochemical  modeling  in¬ 
cludes  an  evaluation  of  ohmic,  activation  and  diffusion 
losses  as  well  as  a  kinetic  model  of  hydrocarbon  reac¬ 
tions,  based  on  most  recent  literature  experiences.  The 
model  has  been  applied  to  a  tubular  geometry,  and  has 
been  calibrated  on  the  available  experimental  data  for  at¬ 
mospheric  and  pressurized  tubular  SOFCs,  showing  the 
capacity  of  predicting  accurately  the  SOFC  operating 
conditions. 

The  model  generates  total  cell  balances  and  internal  cell 
profiles  for  any  relevant  thermodynamic  or  electrochem¬ 
ical  variable,  giving  the  possibility  of  discussing  the  ef¬ 
fects  of  different  operating  conditions  on  the  internal  FC 
behavior. 

A  sensitivity  analysis  applied  to  the  model  parameters  has 
shown  the  effect  of  different  hypothesis  for  the  evaluation  of 
the  cell  losses,  as  well  as  of  the  cell  internal  heat  exchange 
processes  and  of  hydrocarbons  reforming  reactions.  Among 
other  results,  it  is  shown  that  the  importance  of  the  adop¬ 
tion  of  appropriate  parameters  for  the  evaluation  of  activa¬ 
tion  polarization,  as  well  as  the  relevance  of  a  kinetic  model 
for  reforming  reactions,  while  at  the  high  operating  temper¬ 
atures  typical  of  SOFCs,  diffusion  losses  only  slightly  affect 
the  fuel  cell  operation. 

In  future  works,  the  model  will  be  extended  to  other  SOFC 
geometries  and  will  be  applied  to  simulate  SOFC  perfor¬ 
mances  when  integrated  in  more  complex  plant  configura¬ 
tions  with  variable  operating  conditions. 


6  Value  of  pre-exponential  coefficient  and  activation  energy  are  not 
specified  in  Ref.  [17]  and  are  kept  as  by  Table  1. 
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Appendix  A.  Anode  and  cathode  diffusion  polarization 

Reactant  concentration  at  cell  reaction  sites  are  used  in 
Eq.  (20)  for  the  assessment  of  diffusion  polarization.  Calcu¬ 
lation  of  hydrogen,  steam  and  oxygen  molar  fraction  con¬ 
centration  at  cell  reaction  sites  is  carried  out  as  discussed 
below. 

A.L  Anode  reactant  concentration 


Hydrogen  consumption  at  cell  reaction  sites  generates  a 
concentration  gradient  and  a  hydrogen  flux  which  can  be 
estimated  by  the  Fick’s  law  [34]: 


PmD  dXu2 
RvT  dr 


+  Ah2  Aot 


(A.l) 


The  path  of  the  hydrogen  flux  includes  (i)  diffusion  in  the 
gas  mixture  from  the  bulk  flow  composition  (“b”)  to  the  cell 
surface  layer  (“1”)  and  (ii)  diffusion  through  the  cell  porous 
electrode  to  cell  reaction  sites  (“r”). 

The  first  step  is  calculated  as  ordinary  diffusion  in  a  gas 
mixture,  with  a  diffusion  coefficient  [34,35]  equal  to 


1~*h2 

Hi  °H 2,i 


(A.2) 


where  Z)h2,;  is  the  mutual  diffusion  coefficient  of  hydrogen 
in  the  generic  ith  chemical  specie  (/  =  H2O,  CO,  CO2,  CH4, 
N2),  calculated  by  the  Fuller  equation  [24]: 


Dn2,i  = 


0.00143rfL75 


p  m1,/2  rv1/3  +  iV3i2 

'  an ‘v‘\\2  j L'  | |2  T  vs  J 


(A. 3) 


where  Mh2,*  and  v;  can  be  found  in  literature  [24]. 

By  integrating  Eq.  (A.l)  along  the  first  diffusion  step  Lba 
(equal  to  5  mm  according  to  [14])  with  the  further  condi¬ 
tions: 


7tot  =  7r2  +  7h2o  =  0 


i 

IF 


(steam  flux  is  equal  and  opposite  to 
the  hydrogen  flow)  (A.4) 

(A. 5) 


we  have 


iRa  7f 


7FPanZ)H2,m 


dr 


(A. 6) 


which  yields  the  hydrogen  molar  concentration  at  the  elec¬ 
trode  surface: 

/Rg7fLba 


*h2  =  VL  - 


H2  2FPLmD\\2 


(A. 7) 
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The  same  equations  applied  to  steam  diffusion  give 

fRgTfLba 


<0  =  xb 


H2°  +  2FPmDH2o,m 


(A.8) 


The  diffusion  process  inside  the  porous  anode  material  may 
be  calculated  by  the  Knudsen  model  [16]  with  a  diffusion 
coefficient  calculated  as  [34] 

1/2 

+ 


£>h2,p  =  - 


( 


nMH: 


l 


por 


dh2, 


m 


(A. 9) 


where  e,  r,  rpore  are  given  in  Table  1  according  to  [16]. 

Integrating  Eq.  (A.l)  along  the  anode  length  Lan  allows  to 
calculate  hydrogen  and  steam  molar  fraction  at  cell  reaction 
sites: 

iRg  7f  Lan 

h2  ~~  '  - 


XL-  =  xL  - 


X\ho  =  xu2o  + 


H2  2  FPanDR2,p 

iRg  A  Tan 


(A.  10) 


2EPan^H20,p 


(A.  11) 


A.  2.  Cathode  reactant  concentration 


The  same  method  used  for  the  anode  may  be  applied 
to  oxygen  consumption  at  the  cathode:  oxygen  flow  is  ex¬ 
pressed  as 


T_  _  ^catA*  dX()2 

J<h  ~  ~  R„T  dr 


(A.  12) 


The  diffusion  coefficient  of  oxygen  into  air  (assumed  as  an 
O2-N2  mixture)  can  be  estimated  by  the  Fuller  correlation: 


Ao2,n2  = 


0.00 1437^  75 


p  M1^2  \v 4-  v ^312 
^an^o2jN2L^02  "I"  J 


(A.  13) 


N2 


Eq.  (A.  12)  can  be  integrated  with  the  complementary  con¬ 
ditions: 


J«*  ~  J°2  ~  4f 
Lbc  =  r  cat  —  rT 


l 


x 


o2  dX0 


b 

02 


T -L 


Lbc 


iRgTc 


^02  —  A  JO  4E7>cat/)o2,N2 


dr 


Ao2- 

yielding  the  oxygen  fraction  at  the  cathode  surface: 

,  h  /  RaTcLbci 

Xk,  =  l  +  (^-l)exp/  gc 


(A.  14) 
(A.  15) 

(A.  16) 


4FPcat/)o2,N2 


(A.  17) 


Diffusion  through  the  porous  cathode  structure  is  modeled 
with  Knudsen  diffusion  coefficient: 

1/2 

+ 


Ao2,p  —  ^ 


( 


nMa 


1 


4r por  \27?gTs 


Ad2,n2 


(A.  18) 


which  together  with  Eq.  (A.  12),  given  cathode  length  Lca t- 
yields  the  oxygen  fraction  at  cell  reaction  site: 


r  1  (  RgT s^cat*  \ 

XL  =  1  +  (Xk  -  1)  exp  — g  ) 

02  02  \  AFP  cat£>o2ip  / 


(A.  19) 
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